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Abstract 

This  paper  addresses  the  development  of  a  nonlinear  control  design  for  attenuating  structural  vibrations 
using  magnetostrictive  transducers  operating  in  nonlinear  and  highly  hysteretic  operating  regimes.  We  consider 
as  a  prototype  a  thin  plate  subjected  to  exogenous  pressure  waves  and  controlled  via  Terfenol-D  transducers  at 
the  plate  edges;  however  the  methodology  is  sufficiently  general  to  encompass  a  wide  range  of  structures  and 
magnetic  transducer  designs.  Hysteresis  inherent  to  the  transducer  materials  is  quantified  using  a  homogenized 
energy  framework  and  the  resulting  nonlinear  constitutive  relations  are  used  to  construct  a  PDE  representation 
and  corresponding  finite  dimensional  model  of  the  structural  system.  We  employ  optimal  control  theory  to 
construct  nonlinear  open  loop  control  inputs  which  accommodate  the  hysteresis  inherent  to  the  transducers  but 
are  not  robust  with  regard  to  unmodeled  dynamics  or  disturbances.  Robustness  is  incorporated  by  employing 
perturbation  techniques  to  provide  linear  feedback  laws  acting  on  measured  disturbances.  As  illustrated  via 
numerical  examples,  the  resulting  hybrid  control  design  provides  excellent  control  authority  and  robustness  for 
transducers  operating  in  hysteretic  and  nonlinear  regimes. 
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1.  Introduction 

The  active  and  passive  attenuation  of  vibrations  in  aeronautic,  aerospace,  automotive,  and  industrial  systems 

constitutes  a  fundamental  component  of  structural  and  structural  acoustic  design.  These  vibrations  can  be  due 

to  a  wide  range  of  exogenous  inputs  including  adjacent  machinery,  environmental  inputs,  or  impinging  acoustic 

or  fluid  fields,  and  can  be  periodic,  quasiperiodic,  or  random  in  nature.  For  certain  systems,  passive  damping 
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techniques  prove  adequate  and  are  advantageous  due  to  their  simplicity  and  robustness.  However,  passive 
vibration  control  often  comes  at  the  cost  of  added  weight  and  size,  and  this  approach  often  proves  inadequate  in 
moderate  to  high  frequency  regimes.  For  such  applications,  active  or  semi-active  vibration  control  is  required. 
Consideration  of  active  designs  is  further  motivated  by  increased  development  and  utilization  of  highly  flexible 
and  lightweight  structural  components  having  little  inherent  damping  —  e.g.,  composites  and  polymers  -  and 
increasingly  stringent  design  criteria  and  operating  regimes. 

Over  the  past  two  decades,  a  range  of  smart  materials  —  including  piezoelectric  compounds,  relaxor  fer¬ 
roelectric  materials,  magnetostrictive  compounds,  and  shape  memory  alloys  (SMA)  —  have  emerged  as  viable 
alternatives  for  both  passive  and  active  vibration  control.  Their  advantages  arise  from  number  of  factors  in¬ 
cluding  multifunctionality  (e.g.,  both  actuator  and  sensor  capabilities),  large  force  or  strain  generation,  high 
frequency  and  broadband  actuator  capabilities  (excluding  SMA),  and  the  potential  for  minimal  weight  increase. 
At  low  drive  levels,  the  dynamics  of  ferroelectric  and  ferromagnetic  actuators  can  often  be  adequately  approx¬ 
imated  using  linear  constitutive  relations  which  subsequently  form  the  basis  for  linear  structural  models  and 
control  designs.  In  moderate  and  high  drive  regimes,  however,  inherent  hysteresis  and  constitutive  nonlinearities 
sufficiently  dominate  material  behavior  to  require  inclusion  in  models  and  model-based  control  designs. 

As  detailed  in  [40],  there  exists  a  wide  range  of  techniques  for  modeling  the  hysteresis  inherent  to  ferro¬ 
electric  (e.g.,  PZT,  PMN  below  the  glass  transition  temperature),  ferromagnetic  (e.g.,  Terfenol-D,  steel)  and 
ferroelastic  (e.g.,  SMA)  compounds.  Three  modeling  classes  which  provide  unified  characterization  frameworks 
for  ferroelectric,  ferromagnetic,  and  ferroelastic — collectively  termed  ferroic  -  compounds  are  domain  wall  mod¬ 
els  [25,  28,  42],  Preisach  models  [22,  35,  37],  and  homogenized  energy  models  [38,  41,  44,  45].  The  domain  wall 
models  are  efficient  to  implement  but  require  a  priori  knowledge  of  turning  points  to  guarantee  closure  of  biased 
minor  loops.  For  feedback  design,  this  can  prove  problematic  since  turning  points  are  dictated  by  state  mea¬ 
surements  or  estimates  for  the  underlying  structure.  For  certain  operating  regimes,  the  utility  of  a  model-based 
optimal  control  design  utilizing  a  domain  wall  model  is  illustrated  in  [39]  —  for  general  tracking  and  regulation 
application  however,  feedback  control  laws  exploiting  domain  wall  models  should  be  used  with  care  due  to  the 
inherent  nonclosure  property.  Preisach  models  were  originally  developed  for  magnetic  materials  [17,  35]  and 
have  been  exploited  for  model-based  control  design  for  all  three  classes  of  compounds.  The  primary  advantage 
of  this  technique  is  also  one  of  its  primary  limitations,  namely  its  general  mathematical  nature.  For  systems  in 
which  the  physics  is  poorly  understood  or  difficult  to  quantify,  this  generality  is  a  definite  advantage.  However, 
it  is  difficult  to  correlate  parameters  in  the  framework  with  measured  quantities,  and  extensions  to  the  theory 
are  required  to  accommodate  noncongruency,  reversibility,  nonclosure  due  to  relaxation  mechanisms,  and  stress 
or  temperature  dependencies  [15,  16,  40].  In  this  paper,  we  employ  the  homogenized  energy  framework  due  to 
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its  energy  basis  and  flexibility  with  regard  to  numerous  operating  conditions.  Details  regarding  this  framework, 
including  its  relation  to  certain  extended  Preisach  formulations,  can  be  found  in  [40]. 

The  majority  of  control  designs  utilizing  smart  material  transducers  have  focused  on  linear  models  which 
are  effective  in  low  to  moderate  drive  level  regimes  [4,  6,  32,  33,  34].  Within  the  context  of  nonlinear  hysteresis 
models,  there  are  two  primary  approaches:  (i)  linear  control  design  utilizing  model  inverses  employed  as  inverse 
filters  to  approximately  linearize  the  transducer  response,  and  (ii)  nonlinear  control  design.  An  overview  of  the 
former  approach  using  a  piecewise  linear  Preisach  representation  can  be  found  in  [46]  and  the  use  of  this  technique 
for  robust  control  design  utilizing  the  homogenized  energy  framework  is  detailed  in  [29] .  This  technique  has  the 
advantage  that  it  permits  linear  adaptive,  optimal,  classical  or  robust  control  design.  The  primary  disadvantage 
resides  in  the  fact  that  inputs  to  the  model  inverse  often  are  more  difficult  to  interpret  physically  than  direct 
inputs  to  the  transducer.  In  category  (ii),  nonlinear  control  theory  is  used  to  ascertain  inputs  to  the  actuator 
which  yield  the  desired  system  response  —  e.g.,  effective  vibration  attenuation  or  high  accuracy  tracking.  This 
category  includes  the  domain  wall-based  optimal  design  in  [39]  as  well  as  nonlinear  designs  of  the  type  detailed 
in  [21,  3,  23].  Schematics  illustrating  the  two  approaches  are  provided  in  Figure  1. 


In  this  paper,  we  develop  a  nonlinear  optimal  control  technique  for  magnetostrictive  transducers  which 
utilizes  the  homogenized  energy  framework  to  quantify  the  hysteretic  map  between  input  magnetic  fields  H 
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Figure  1:  (i)  Linear  control  design  employing  an  inverse  filter,  and  (ii)  nonlinear  control  design. 
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and  the  magnetization  M  and  quadratic  magnetization  -  strain  (M  —  e)  relation.  Magnetostrictive  actuators 
have  been  considered  for  a  number  of  vibration  attenuation  applications  including  the  reduction  of  chatter  in 
machining  operations,  vibration  attenuation  in  transmissions  and  axle  differentials,  and  damping  of  large  flexible 
structures  [12,  13,  36].  To  provide  a  prototype  for  which  the  structural  physics  is  well-understood  but  finite¬ 
dimensional  models  are  relatively  high-order,  we  consider  the  control  of  vibrations  in  a  cantilever  plate  using 
edge-mounted  Terfenol-D  transducers  as  depicted  in  Figure  2.  This  simplified  construct  permits  us  to  focus  the 
discussion  on  attributes  of  the  nonlinear  control  design  for  a  2-D  structure  without  the  more  complicated  system 
models  required  for  the  previously  mentioned  applications.  In  this  context,  we  illustrate  the  effectiveness  of  the 
technique  for  systems  having  over  150  degrees  of  freedom  which  provides  an  indication  that  the  control  method 
will  be  feasible  when  applied  to  finite  element  models  for  more  complex  applications  and  designs. 

The  hybrid  nonlinear  control  design  is  comprised  of  an  open  loop  component,  derived  from  nonlinear,  finite¬ 
dimensional  control  theory,  and  a  feedback  component  obtained  through  perturbation  techniques.  The  open  loop 
control  accommodates  the  hysteresis  and  constitutive  nonlinearities  inherent  to  magnetostrictive  transducers 
but  is  not  robust  with  regard  to  unmodeled  dynamics  or  uncertainties  in  operating  conditions.  Robustness  is 
obtained  by  linearizing  the  system  about  the  optimal  trajectory  and  control,  and  feeding  back  on  perturbations 
from  the  open  loop  system  —  e.g.,  see  [5]  for  analysis  of  perturbation  control  for  general  nonlinear  systems  or  [27] 
for  certain  algorithms.  This  approach  has  the  advantage  that  the  resulting  perturbation  system  is  linear  with 
quadratic  constraints  and  hence  LQR  theory  can  be  applied  to  facilitate  implementation.  The  implementation 
procedure  can  be  summarized  as  follows,  (i)  The  nonlinear  two-point  boundary  value  problem  required  to 
obtain  an  open  loop  control  is  solved  off-line  and  the  resulting  optimality  system  is  stored,  (ii)  During  numerical 
or  experimental  implementation  of  the  method,  perturbations  from  the  optimality  system  are  measured  and 
feedback  via  either  LQR  or  classic  methods  to  compensate  for  unmodeled  dynamics  or  uncertainties  in  operating 
conditions.  As  demonstrated  through  numerical  examples,  the  resulting  hybrid  control  method  is  robust  and 
efficient  to  implement.  Hence  it  presents  an  efficient  alternative  to  dynamic  programming  solutions  [5]  or 
state-dependent  Ricatti  solutions  [1]  for  real-time  implementation  of  nonlinear  control  laws. 

The  results  presented  here  differ  from  those  of  [39]  in  two  fundamental  aspects:  (i)  the  homogenized  energy 
framework  provides  significantly  more  flexibility  for  control  design  using  a  range  of  smart  material  actuators 
than  the  domain  wall  model  employed  in  [39],  and  (ii)  the  development  and  illustration  of  the  theory  in  the 
context  of  a  2-D  plate  model  rather  than  the  1-D  beam  model  employed  in  [39]  demonstrates  the  potential  for 
implementing  the  method  with  larger  scale  FE  models  used  to  model  the  dynamics  of  more  general  vibration 
attenuation  systems. 

The  development  of  the  optimal  control  method  is  presented  as  follows.  Section  2  summarizes  magne- 
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tostrictive  material  behavior  and  the  homogenized  energy  model.  The  magnetically  actuated  plate  structure  is 
developed  in  Section  3.  In  Section  4,  the  optimal  control  problem  is  discussed.  Linear  optimal  control  results 
illustrate  the  need  for  nonlinear  control  when  field  levels  reach  moderate  to  high  regimes.  The  nonlinear  optimal 
control  method  is  then  developed  and  validated  numerically  for  the  cases  of  free  and  forced  vibration.  Finally 
perturbation  control  is  implemented  to  improve  robustness  to  operating  uncertainties. 

2.  Material  Behavior  and  Hysteresis  Model 

Magnetostrictive  materials  are  a  special  class  of  magnetic  compounds  that  undergo  a  shape  change  when 
exposed  to  a  magnetic  field.  The  deformation  is  a  manifestation  of  local  magnetic  moments  that  align  with 
the  applied  field.  Under  small  fields,  regions  of  like  magnetization  called  domains  move  in  an  approximately 
reversible  manner  —  creating  small  strains  and  minimal  hysteresis.  As  fields  increase,  irreversible  domain  wall 
motion  and  rotation  of  domains  develop  creating  significant  hysteresis  and  nonlinearity. 

The  pseudobinary  compound  Terfenol-D  (Tbo.3Dyo.7Fe1. 9)  is  increasingly  employed  in  actuator  design  due 
to  its  large  force  and  strain  response  under  a  magnetic  field  [6,  9].  Magnetostriction  in  the  approximately  linear 
range  is  on  the  order  of  500  microstrain,  whereas  strain  in  excess  of  1000  microstrain  can  be  achieved  in  the 
nonlinear  range.  Certain  transducer  designs  can  generate  forces  up  to  550  N  with  broadband  capability  from 
DC  up  to  20  kHz  [20], 

The  magnetostrictive  transducer  employed  in  the  structural  control  problem  is  described  in  detail  by  Dapino 
et  al.  [13].  As  shown  in  Figure  2(a),  the  actuator  primarily  consists  of  a  Terfenol-D  rod,  a  surrounding 
wire  solenoid,  a  permanent  magnet,  and  a  spring  washer /compression  bolt  assembly  to  prestress  the  rod.  A 
current  is  applied  to  the  solenoid  which  generates  a  magnetic  field.  The  magnetic  field  produces  magnetic  flux, 
magnetization,  and  strain  in  the  Terfenol-D  rod.  The  rod  is  prestressed  to  reduce  tensile  loads  during  operation 
and  to  increase  actuator  displacement  [18].  The  permanent  magnet  serves  to  bias  the  magnetization  on  the 
Terfenol-D  rod  such  that  bi-directional  strains  are  produced  by  a  time-varying  magnetic  field  with  zero  bias 
and  to  provide  uniform  flux  patterns  [11].  Bi-directional  strain  can  also  be  achieved  by  biasing  the  magnetic 
field,  but  this  results  in  substantial  power  losses  from  ohmic  heating  which  necessitates  use  of  liquid-cooled 
transducers.  Trade-offs  between  size  and  weight  of  permanent  magnets  and  power  loss  from  ohmic  heating  must 
be  considered  for  a  given  application. 

2.1  Homogenized  Energy  Model 

A  one-dimensional  magnetostrictive  constitutive  law  is  used  to  model  the  Terfenol  transducers.  The  homog¬ 
enized  energy  model  focuses  on  the  nonlinear,  hysteretic  H  —  M  behavior  while  assuming  linear  stress-strain 
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Figure  2:  (a)  Terfenol-D  transducer  design,  and  (b)  transducer  configuration  for  attenuating  transverse  plate  vibrations. 


behavior.  A  brief  summary  detailing  the  constitutive  model  is  presented  here.  Details  regarding  the  model 
development  are  found  in  [40,  41]. 

To  summarize,  a  Gibbs  relation  at  the  mesoscopic  length  scale  is  used  to  quantify  the  local  magnetization  M 
as  a  function  of  the  applied  field  H  for  regimes  in  which  relaxation  mechanisms  may  be  negligible  or  significant. 
The  model  is  extended  to  macroscopic  scales  through  stochastic  homogenization  techniques.  The  Gibbs  relation 
is 


G{M)  =  T(M)  -  no HM 

where  T ( M ]  is  the  Helmholtz  energy  [40]  and  /io  is  the  permeability  of  free  space. 
As  detailed  in  Chapter  7  of  [40],  the  constitutive  relations  are 


(1) 


a(t)  =  EMe(t )  -  ar  ([M(H)}  ( t )  -  M0)2 

(2) 

f-OC  pQO  v  ' 

\M(H)]  (t)=  /  Vl{Hc)v2{H!)  [M  (H  +  Hr,  Hc ,  0]  (t)dHIdHc 

J  0  J — oo 

where  a(t)  is  the  uniaxial  stress,  e{t)  is  the  elastic  strain,  EM  is  the  elastic  modulus  at  constant  magnetization,  ai 
is  the  magnetostrictive  coefficient,  M  is  the  magnetization  and  Mq  includes  magnetization  from  the  permanent 
magnet  and  the  initial  internal  state  of  the  material.  Here,  Hi  is  the  interaction  field,  Hc  is  the  coercive  field 
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and  £  denotes  the  initial  distribution  of  magnetic  variants.  The  stochastic  homogenization  technique  used  to 
formulate  the  macroscopic  magnetization  M  in  terms  of  local  values  M  is  based  on  the  assumption  that  the  local 
coercive  field  Hc  and  interaction  field  Hi  are  manifestations  of  underlying  distributions  rather  than  constants  to 
accommodate  a  myriad  of  local  material  inhomogeneities  such  as  intergranular  residual  stress,  impurities,  and 
certain  anisotropies  .  The  corresponding  densities  are  designated  by  Vi(Hc)  and  u2(Hi).  As  detailed  in  [41], 
one  choice  for  these  densities  is 


MHMHj)  =  creae-N^J/HV*?^  (3) 

where  Hc  is  the  average  coercive  field,  c  quantifies  the  coercive  field  variability,  b  is  the  variance  of  the  interaction 
field,  and  ci  and  c2  are  scaling  parameters.  The  proposed  densities  are  implemented  to  reduce  parameter 
estimation.  Model  comparison  to  experimental  results  can  be  found  in  [41]. 

For  many  operating  regimes,  thermal  activation  mechanisms  such  as  magnetic  after-effects  and  accommoda¬ 
tion  [24]  significantly  affect  the  constitutive  behavior.  These  effects  are  incorporated  in  the  framework  through 
Boltzmann’s  relations 


H(G)  =  Ce~GV/kT  (4) 

which  balances  the  Gibbs  and  relative  thermal  energies.  Here  G  is  the  Gibbs  energy,  V  is  a  representative  volume 
element,  k  is  Bolztmann’s  constant,  and  T  is  temperature.  The  constant  C  is  specified  to  ensure  integration  to 
unity. 

In  this  case,  the  local  magnetization  is  defined  by 


M  =  x+(M+)  +  x  _(M_>  (5) 

where  x+  and  X-  respectively  denote  the  volume  fraction  of  moments  having  positive  and  negative  orientations. 
The  variables  (M+)  and  (M_)  are  the  average  magnetizations  corresponding  to  the  volume  fractions  of  moments. 
The  differiential  equations  governing  evolution  of  x+  and  X-  is  based  on  material-dependent  parameters  that 
define  the  likelihoods  that  moments  switch  from  positive  to  negative,  and  conversely.  These  likelihoods  account 
for  the  observed  thermal  relaxation  mechanisms.  Details  describing  the  governing  equations  are  given  in  [41]. 

In  certain  operating  regimes,  thermal  relaxation  is  negligible  and  in  which  case  the  relation  (5)  limits  to  the 
piecewise  linear  relation 


M(H  +  Hr,  Hc.  0  =Xm(H  +  Hi)  +  MRS(H  +  Hr  Hc,£) 


(6) 
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where  Mu  is  the  magnitude  of  the  local  magnetization  variant.  The  variable  5  is  equal  to  one  if  the  magneti¬ 
zation  variant  is  in  the  positive  direction  and  negative  one  if  oriented  in  the  negative  direction.  The  magnetic 
susceptibility  is  defined  by  Xm- 

The  local  magnetization  M f>  given  in  (6)  will  switch  when  magnetic  variants  diametrically  opposed  to  the 
effective  field  ( He  =  H  +  Hj)  reach  the  coercive  field.  The  switching  behavior  is  modeled  by  introducing  a 
semi-infinite  set  of  magnetic  variants  that  correspond  to  the  distribution  of  effective  fields  and  coercive  fields. 
Details  regarding  the  numerical  implementation  of  the  switching  behavior  and  the  numerical  integration  of  (2) 
are  provided  in  [40,  45].  The  constitutive  response  predicted  by  the  homogenized  free  energy  model  is  illustrated 
in  Figure  3.  The  local  magnetization  given  by  (6)  has  been  used.  Although  this  is  a  linear  piecewise  function, 
the  typical  nonlinearities  associated  with  magnetostrictive  materials  are  captured  by  the  local  relation  for  M. 

For  purposes  of  developing  the  structural  model  and  control  design,  the  local  magnetization  given  by  (6)  is 
implemented  where  thermal  activation  mechanisms  are  neglected. 


Magnetic  field  H  (kA/m) 

(b) 


Figure  3:  (a)  Macroscopic  magnetization  M  versus  magnetic  field  H  computed  using  (2),  and  (b)  microstrain  /re  versus 
magnetic  field  H  computed  from  (2)  under  zero  external  stress. 


3.  Structural  Model 


The  structure  under  consideration  consists  of  Terfenol  transducers  positioned  along  two  fixed  edges  of  a 
thin  plate  as  depicted  in  Figure  2.  Eight  transducers  are  oriented  in  four  pairs  with  each  pair  attached  to  the 
plate  by  a  moment  arm  that  is  perpendicular  to  the  plate.  Diametrically  out-of-phase  currents  applied  to  each 
actuator  couple  generate  a  moment  on  the  plate. 

3.1  Plate  With  Nonlinear  Actuators 

The  thin  plate  structure  is  modeled  using  classical  plate  theory  [47].  The  equation  of  motion  describing 
transverse  plate  displacement  w  is  decoupled  from  in-plane  displacements  by  focusing  on  small  strain  induced 
by  plate  vibration.  Since  the  control  problem  is  concerned  with  transverse  plate  vibration,  we  focus  only  on  the 
dynamic  equation  for  w. 

The  balance  of  forces  and  moments  for  the  plate  result  in  the  strong  form  of  the  equation  of  motion, 

d2w  d2Mint  d2M d2Mint  d2M™ag  d2M™ag 
Pv~dV  ~d. 2~d^dy  ~  dx2  +  dy 2  +9 


where  pv  denotes  the  density  of  the  plate,  A4™4  and  A4™*  include  the  internal  elastic  and  damping  moment 
components  in  the  x  and  y  directions,  and  A4™*  is  the  in-plane  twisting  moment.  The  external  moments 
generated  by  the  Terfenol  transducers  are  represented  by  A4”los  and  A4™°s  in  the  x  and  y  directions,  respectively. 
No  twisting  moments  are  applied  by  the  transducers,  so  the  component  is  zero.  The  external  disturbance 

load  distributed  along  the  plate  surface  is  denoted  by  g.  A  detailed  description  of  the  equations  describing  the 
moment-displacement  relations  are  given  in  the  Appendix.  Spatial  variations  in  density,  moment  of  inertia  and 
compliance  near  the  Terfenol  transducers  and  connecting  moment  arms  are  assumed  to  be  negligible. 

The  plate  geometry  is  represented  by  the  length  £x,  width  £y  and  thickness  h.  Values  used  in  the  model  are 
given  in  Table  1.  The  plate  is  assumed  to  be  clamped  along  the  fixed  edges  [x,  0]  and  [0,  y\  and  free  along  [lx,  y\ 
and  [x,ty\  as  shown  in  Figure  2(b).  The  cross-sectional  area  of  the  Terfenol  rod  is  denoted  by  Amag  and  we  let 
lr  deonote  the  length  of  the  moment  arm  from  the  actuator  to  the  surface  of  the  plate.  The  location  of  the  four 
connecting  rods  are  designated  by  fij,  i  =  1, . . . ,  4.  Let  -Hi(j)  and  *  =  1, . . , ,  4  respectively  denote  the  field 
inputs  to  the  top  and  bottom  actuators  in  each  pair  and  let  H(t)  =  [i?1(1)(£),  H2(i){t), . . . 
denote  the  complete  set  of  inputs.  The  external  moments  generated  by  the  Terfenol  actuators  can  then  be 
expressed  as 
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(8) 


Mrgm)^,y)  =  (#«(*))  -  M22(i)  (H{i)(t))}Xx{i)(x,y) 

i= 1 

4 

M™a9(H(t),x,  y)  =  -KM  (H(i)  W)  -  Mki)  (tf(i)W)]xy(i)(*.  y)  (9) 

i= 1 

where  ICM  =  a\Amag  (h/2  +  £r)  and  and  M2(i)(t)  are  the  magnetizations  modeled  by  (2).  The  charac¬ 

teristic  functions 


Xx(i){x,y) 


Xv(i)(x,y) 


1  ,  (x,y)  G  Qi,  i  =  1,2 

0  ,  otherwise 

1  ,  (x,y)  G  Qi,i  =  3,4 

0  ,  otherwise 


(10) 


(11) 


designate  that  actuators  1  and  2  along  [0 ,  y\  generate  moments  in  the  x-direction  whereas  actuators  3  and  4 
along  [x,  0]  generate  moments  in  the  y-direction. 

The  control  design  is  simplified  by  reducing  the  number  of  control  inputs  from  eight  to  four.  The  time 
varying  control  inputs  are  computed  for  each  of  the  four  top  transducers  and  the  control  inputs  on  the  bottom 
four  actuators  are  defined  by  The  magnetization  is  then  defined  by  the  relation  Mpj)  = 

M(j)  +  Mo  and  M2p)  =  +  Mq  which  assumes  negligible  differences  between  the  magnetization  driven  by 

the  fields  applied  on  the  top  and  bottom  actuators.  This  yields  the  relation 


4 

Mrgm),x,y)  =  ^M(i)(ff(0(t))xx(o(x,y)  (12) 

2=1 

4 

M™aa(H(t),x,y)  ^Mw(M(i)(t))X!/W(x,y)  (13) 

2=1 

where 


=  -/CM[4M0M(i)  (H{i)(t))}  (14) 

which  is  an  approximation  of  actual  magnetostrictive  material  behavior  when  hysteresis  and  nonlinearities  are 
present.  The  control  law  can  be  improved  by  including  all  eight  transducers  as  individual  inputs,  but  for  the 
purposes  of  demonstrating  the  control  design  the  simpler  four  input  model  is  used. 

As  detailed  in  [40],  the  strong  model  formulation  (7)  can  be  rewritten  in  the  weak  form 


10 


-  Mint^~  -  2Mvnt  °  ^ 

'  dt 2  *  ax2  dccdy 


■  ,a2$\  /■  (  a2$  a2$  \ 

•^r^J  ^  =  Jn  [Mr3g^+M™°W+9<t>)  du  (15) 


where  $  £  Hg(x,y)  and  O  =  [0,4;]  x  [0, 4y]  denotes  the  plate  region. 

3.2  Approximation  Method 

The  plate  displacements  can  be  approximated  using  either  cubic  B-splines  or  cubic  Hermite  elements.  Cubic 
B-splines  are  chosen  for  the  control  problem  primarily  because  cubic  Hermite  elements  require  solving  for 
approximately  twice  the  number  of  unknown  coefficients  which  may  affect  computation  speed  for  real-time 
control.  Details  describing  the  attributes  of  cubic  B-splines  as  well  as  comparisons  to  cubic  Hermite  elements 
are  given  in  Chapter  8  of  [40]. 

The  cubic  B-splines  are  defined  over  the  plate  geometry,  H  =  [0,  £XJ  x  [0 ,£y\.  The  plate  is  discretized  using 
the  points  xm  =  mhx  and  ym  =  nhy  with  hx  =  hy  =  jf-  and  m  =  0, Nx  and  n  =  0,  ...,Ny.  The  cubic 
spline  product  space  is  defined  by 


${x,y)  =  <t>(x)<i>(y)-  (16) 

The  approximate  solution  to  (15)  is  subsequently  given  by 

Nw 

w(t,  x,y)  =  ^2  wk(t)$k(x,  V )  (17) 

k= 1 

where  Nw  =  ( Nx  +  1  )(Ny  +  1)  and  Wk(t)  are  the  coefficients  to  be  determined  through  solution  of  the  weak 
model  formulation. 

Approximation  of  the  infinite  dimensional  model  (15)  yields  the  ODE  system 

Mw  +  Cw  +  Kw  =  bM  (H(f))  +  g  (18) 

where  M,  C,  and  IK  are  the  mass,  damping  and  stiffness  matrices.  The  vector  w  =  [■uq(t),  ...,wnw  (t)]T  represents 
the  coefficients  in  (17)  and  g  consolidates  disturbance  loads.  Details  describing  the  components  of  these  matrices 
and  vectors  are  given  in  the  Appendix  and  in  [40].  As  indicated  in  (54)  in  the  Appendix,  the  (Nw  +  1)  x  4 
matrix  b  specifies  the  spatial  components  of  the  inputs.  The  4x1  input  vector 

M  (H(f))  =  [jM(1)  (H(1)(t))  , . . . ,  A4(4)  (i7(4)(f))]T  (19) 
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designates  the  time-dependent  influence  of  the  field  inputs. 

The  second-order  system  can  be  formulated  as  a  set  of  first-order  equations  for  implementation  in  the  control 
design 


y{t)  =  Ay(t)  +  [-B(u)](i)  +  G{t) 

2/(0)  =  2/o 

where  y(t)  =  [w\ (t) ,  •  •  •  ,  (t),  wi(t),  ■  ■  ■  ,W]yw(t)}.  The  system  matrix  and  input  vector  are  given  by 


(20) 


- 1 

0 

1=1 

_ 1 

0 

,  B(u)  = 

— M-1C 

M-1b 

M(u) 


(21) 


with  a  similar  definition  for  the  disturbance  load  G(t).  The  identity  matrix,  with  dimension  Nw  x  Nw,  is  denoted 
by  I.  The  control  input  u  is  the  magnetic  field  H  applied  to  each  transducer. 


3.3  Structural  Parameters 

Physical  parameters  employed  in  the  control  design  are  summarized  in  Table  1.  The  Terfenol  material 
parameters  implemented  in  the  homogenized  free  energy  model  are  within  the  range  obtained  for  model  fits  to 
an  experimental  transducer  [8].  The  plate  modulus  Ep,  Poisson  ratio  up  and  density  pp  are  typical  for  aluminum. 
The  internal  damping  is  assumed  to  be  linear  and  is  represented  by  cp  while  viscous  air  damping  is  given  by 
7.  The  center  location  of  the  four  actuator  pairs  was  (0.105, 0.2),  (0.105,  0.3),  (0.15, 0.105)  and  (0.25,0.105)  in 
millimeters  where  the  coordinate  origin  was  located  at  the  top  left  corner  of  the  plate  as  shown  in  Figure  2(b). 
The  moment  arm  is  assumed  to  be  rigid  with  cross-sectional  area  denoted  by  Ar.  The  cross-section  area  is 
used  to  determine  the  region  over  which  the  moment  is  applied.  Plate  dynamics  were  sufficiently  resolved  by 
choosing  Nx  =  Ny  =  4  in  the  frequency  range  considered.  The  dimension  of  the  state  vector  y  was  then  50  x  1 
due  to  the  inclusion  of  both  displacement  and  velocity  components.  When  determining  the  convergence  of  the 
numerical  approximations,  however,  we  have  both  simulated  and  run  control  problems  with  state  vectors  up  to 
Nx  =  Ny  =  8  for  a  total  of  162  degrees  of  freedom. 


12 


Tabic  1:  Parameters  for  the  plate  and  Terfenol  transducer. 


Plate 

Terfenol  Transducer 

Ep  =  4.1  x  1010  N/m2 

Em  =  7.00  x  1010  N/m 2 

vp  =  0.345 

ci  =  C2  =  6.1  x  10-5  m/A 

pp  =  2700  kg/m3 

Hc  =  3.3  x  103  A/m 

cp  =  2.5  x  105  Ns/m2 

c=  0.4 

7  =  0.18  Ns/m2 

b  =  1.5  x  104  A/m 

ix  =  0.4  m 

M0  =  6.618  x  104  A/m 

£y  =  0.6  m 

ar  =  0.006  N/A2 

h  =  0.0016  m 

Amag  =  0.0064  m2 

tr  =  2.54  cm 

Xm  =  1 

Ar  =  1  cm2 

4.  Control  Design 

The  optimal  control  problem  is  first  summarized  to  elucidate  the  technique  used  in  developing  the  nonlinear 
control  design  [5,  26,  27,  31].  The  following  performance  index 

J(u)  =  \yT(tf)^fy(tf)  +  £  [ H(y,u,t )  -  A T(t)y\  dt  (22) 

is  used  to  develop  the  optimal  control  design,  where  the  positive  definite  matrix  11/  penalizes  large  terminal 
values  of  the  state,  H(y,u,t)  is  the  Hamiltonian,  and  A (t)  €  R2 ,V”J  is  a  set  of  Lagrange  multipliers. 

The  Hamiltonian  is 


H{y,  A,  u,  t)  =  L(y,  u,  t)  +  XT  [Ay(t)  +  [B(u)](t)  +  G(t)] 


=  \  [yT(t)Qy(t)  +  uT{t)Ru(t )]  +  \T  [Ay {t)  +  [B(u)]{t)  +  G(t )] 


(23) 


where  the  Lagrangian  L  includes  the  penalties  on  the  states  and  inputs  through  the  semi-definite  maxtrix  Q 
and  the  positive  definite  matrix  R. 

The  optimal  control  problem  requires  solution  of  the  two-point  boundary  value  problem  governed  by  (20) 
and  the  adjoint  or  Lagrange  multiplier  condition 


dX  _  dH 
dt  dy 

Htf)  =  n/!/0/)- 


(24) 
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The  optimal  control  is  determined  from  the  stationary  condition, 


dH 

du 


=  0 


which  results  in  the  nonlinear  optimal  control  input 


u*(t)  = 


(25) 


(26) 


4.1  Control  Parameters 

The  penalties  on  the  states  and  inputs  are  adjusted  in  the  linear  and  nonlinear  control  simulations  to  assess 
the  effect  of  the  hysteretic  constitutive  behavior  of  the  magnetostrictive  transducers  on  vibration  attenuation. 
Values  for  each  case  are  given  in  Table  1.  In  the  linear  simulations,  two  sets  of  values  are  used  to  show 
control  performance  when  nonlinearity  and  hysteresis  are  negligible  or  significant.  The  values  used  in  the  linear 
simulations  as  well  as  the  open  loop  nonlinear  control  and  perturbation  control  are  given  in  Table  2.  The  matrix 
Q  in  (23)  is  given  by 


Q  = 


giK  0 
0  <72  M 


The  penalties  on  the  inputs  are  given  by  Rij  =  where  S  is  the  Kronecker  delta  symbol  with  i,j  =  1  to  4 
for  the  inputs  on  the  four  Terfenol-D  transducer  couples.  The  penalty  on  the  final  state  is  defined  by  11/  =  (73II 
where  II  is  the  solution  to  the  algebraic  Ricatti  equation  described  in  the  following  section. 

Table  2:  Parameters  used  to  penalize  the  states  and  inputs  in  the  linear  and  nonlinear  control  simulations.  The  penalty 

on  the  final  state  is  53  =  5  x  10-8. 


Linear  Control 

Linear  Constitutive  Behavior 

Linear  Control 

Nonlinear  Constitutive  Behavior 

Nonlinear  Control 

Perturbation  Control 

q4  =  5.00  x  101 

<?i  =  5.00  x  101 

qi  =  5.00  x  102 

<7i  =  5.00  x  104 

q2  =  5.00  x  10° 

q2  =  5.00  x  10° 

q2  =  5.00  x  101 

q2  =  5.00  x  103 

n  =  5.00  x  10"6 

n  =  1.25  x  10~6 

n  =  2.50  x  10“6 

n  =  2.50  x  10-6 

r2  =  3.33  x  10"6 

r2  =  8.33  x  10-7 

r2  =  1.67  x  10~6 

r2  =  1.67  x  10"6 

r3  =  2.50  x  10"6 

r3  =  6.25  x  10-7 

r3  =  1.25  x  10~6 

r3  =  1.25  x  10"6 

r4  =  5.00  x  10"6 

r4  =  1.25  x  10~6 

r4  =  2.50  x  10~6 

r4  =  2.50  x  10"6 
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4.2  Linear  Optimal  Control 


Linear  theory  provides  acceptable  control  when  input  currents  to  the  solenoid  are  small.  It  has  been  experi¬ 
mentally  shown  that  a  nearly  linear  relation  exists  between  magnetic  held  and  strain  in  the  Terfenol  transducer 
when  the  held  is  small.  In  this  situation,  an  approximate  model  can  be  attained  through  linearization  about 
some  biased  magnetization.  In  the  linear  model,  the  biased  magnetization  Mq  includes  the  fully  magnetized 
internal  state  of  the  material  and  the  magnetization  from  the  permanent  magnet.  This  results  in  a  linear  input 
operator  where  the  linear  H  —  M  relation  =  xmH(i)(t)  is  been  substituted  into  (14)  where  the  macro¬ 

scopic  magnetic  susceptibility  is  xm  =  0.72.  This  value  is  determined  numerically  from  (2)  for  small  held  inputs. 
In  this  case,  the  moment  from  each  actuator  couple  is 


=  -/CM[4M0xmFw(t)] 

which  is  substituted  into  (12)  and  (13)  to  obtain  the  total  applied  moment. 
Due  to  the  linearity  of  the  constitutive  law  used  here,  the  input  operator  is 


where  the  4x1  vector 


B  =  —AK,M  MoXm 


0 

M-1b 


H  (t) 


(27) 


(28) 


H(i)=  [%)(*),...,  H(4)(t)]T 

represents  the  magnetic  held  on  each  actuator  couple. 

With  this  approximation,  the  corresponding  hrst-order  system  is 

y(t)  =  Ay(t)  +  Bu(t)  +  G(t) 

2/(0)  =  2/o  • 

The  state  constraint  in  (30)  and  adjoint  condition  in  (24)  yields  the  optimality  system 


y{t) 

A 

~BR-1Bt 

y(t ) 

G(t) 

= 

+ 

i 

HO 

-Q 

-AT 

i 

c-+- 

1 _ 

0 

(29) 


(30) 


(31) 


2/(*o)  =  2/o 

x(tf)  =  n/2/(t/)  • 
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Computing  the  optimal  control  requires  solution  of  the  two-point  boundary  value  problem  defined  by  the 
optimality  system.  In  the  linear  case,  a  fundamental  solution  matrix  can  be  obtained  by  solving  the  differential 
Ricatti  equation 


-fi  =  ATn  +  ha  -  ubr-1btu  +  q 

n(t/)  =  iif. 

The  optimal  control  input  is  then  defined  by, 


(32) 


u*(t)  =  -R  1  BT[U(t)y(t)  -  r(f)] 


where  the  variable  r(t)  €  IR2Ar' 


is  the  solution  to  the  auxiliary  differential  equation 


(33) 


r(t)  =  -  [A  -  BR~1BtU]T  r(t)  +  ILG(t) 
r(tf)  =  0. 


(34) 


The  linear  control  theory  can  be  further  simplified  through  the  assumption  that  the  disturbance  force  G  is 
periodic  on  the  time  interval  [toAf]-  In  this  case,  the  appropriate  performance  index  is  given  by 


J(u)  =  [  \  [ yT(t)Qy(t )  +  uT(t)Ru(t)]  dt. 
Jt0  1 


(35) 


This  leads  to  a  steady-state  controller  gain  where  the  fundamental  solution  matrix  is  now  governed  by  the 
algebraic  Ricatti  equation  where  II  =  0  in  (32).  The  resulting  control  input  is 


u*(t)  =  - R~1BT[Ily(t )  -  r(t)].  (36) 

Details  regarding  this  approach  are  given  in  [2,  14]. 

The  control  trajectory  determined  by  (36)  is  used  in  the  following  numerical  examples  to  illustrate  the  need 
for  nonlinear  control  when  the  input  field  reaches  a  magnitude  that  induces  hysteresis. 

4.2.1  Numerical  Example  —  Truncated  External  Force 

We  first  consider  the  case  in  which  the  plate  is  excited  with  an  external  force  for  a  short  time  interval  and 
then  the  force  is  subsequentially  set  to  zero.  The  disturbance  load  G(t)  was  spatially  discretized  and  applied 
uniformly  to  the  plate  with  magnitude  25  N  by  a  summation  of  five  equally  weighted  frequencies  (/e  =10,  12, 
19,  26,  37  Hz)  as  given  by 
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25Ef=i[sm(2 Tr/fi)]  ,  t  <  0.45 


(37) 


g{t,x,y ) 


0  , 


0.45  <  t  <  2.5 


where  E)"=i  [sin(27r/®t)]  =  [sin(27rl0f)  +  •  •  •  +  sin(27r37f)]. 

The  uncontrolled  and  controlled  displacements  at  the  plate  tip  [£x,£y\  are  plotted  in  Figure  4(a).  The 
penalties  on  the  states  Q  and  inputs  R  are  adjusted  to  ensure  the  constitutive  response  computed  using  (2)  is 
linear.  Marginal  vibration  control  is  obtained  under  small  input  fields.  The  linear  H  —  M  response  corresponding 
to  the  input  control  for  each  actuator  pair  is  illustrated  in  Figure  4(b). 

In  Figure  5,  the  linear  control  is  increased  by  adjusting  Q  and  R  (see  Table  2)  to  illustrate  the  effect  of 
hysteresis  on  controlling  vibration.  The  increase  in  magnetic  field  creates  hysteresis  which  is  not  accounted 
for  in  the  linear  control  design.  When  this  occurs,  the  vibration  control  initially  improves,  but  subsequently 
degrades  due  to  the  phase  shift  induced  by  the  hysteresis  illustrated  in  Figure  5(b). 


4.2.2  Numerical  Example  Periodic  External  Disturbance  Force 

Secondly,  we  consider  the  case  in  which  the  plate  is  continuously  excited  with  the  external  disturbance 
force  (37)  (as  shown  in  Figure  2(b))  during  the  entire  simulation.  The  control  is  again  applied  at  fo=0.45  s  and 


(a) 


(b) 


Figure  4:  Performance  of  the  linear  feedback  control  law  using  the  controller  gain  in  (33).  The  resulting  plate  displacement 

for  the  uncontrolled  ( - )  and  controlled  response  (^— )  are  shown  in  (a).  The  relationship  between  magnetic  field  and 

magnetization  is  shown  in  (b)  for  each  actuator  pair. 


17 


Figure  5:  Linear  feedback  control  law  at  levels  where  the  transducer  model  exhibits  hysteresis.  In  this  simulation,  Q  and 

R  were  adjusted  to  increase  the  input  field.  The  resulting  plate  displacement  for  the  uncontrolled  ( - )  and  controlled 

response  (^— )  are  shown  in  (a)  and  the  corresponding  hysteretic  M-H  behavior  is  shown  in  (b). 


vibration  attenuation  is  marginal  as  shown  in  Figure  6.  Additionally,  the  disturbance  forces  induce  nonlinearities 
in  the  magnetostrictive  actuators  from  larger  field  inputs  necessary  to  minimize  vibration  which  is  not  properly 
accounted  for  in  the  linear  control  design. 

4.3  Nonlinear  Control  Method 

The  nonlinear  control  method  previously  discussed  by  Smith  [39]  is  applied  to  the  Terfenol  plate  structure. 
Key  equations  describing  the  numerical  technique  are  summarized  to  elucidate  the  approach  used.  The  input 
operator  B(u)  contains  nonlinearities  from  the  magnetostrictive  constitutive  law  which  does  not  lead  to  a 
fundamental  matrix  solution;  therefore  an  efficient  solution  in  terms  of  the  Riccati  equation  is  not  possible. 
This  issue  is  resolved  by  solving  the  optimality  system  numerically  while  ensuring  the  boundary  conditions  at 
the  initial  and  final  times  are  satisfied.  The  first  order  system  of  equations  is  formulated  as, 

z(t)  =  F(t,z)  (38) 

where  z  =  [y,  A]T  and 
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Figure  6:  Performance  of  the  linear  feedback  control  law  in  the  presence  of  a  periodic  external  force  of  25  N.  The  resulting 

plate  displacement  for  the  uncontrolled  ( - )  and  controlled  response  (^— )  are  shown  in  (a).  The  relationship  between 

magnetic  field  and  magnetization  is  shown  in  (b)  for  each  actuator  pair. 


F(t,z) 


Ay(t)  +  [-B(u)](t)  +  G(t) 
-AT\(t)  -  Qy(t ) 


(39) 


y(to)  =  j/o  ,  a (tf)  =  n fy(tf) . 


The  system  of  equations  given  by  (38)  are  solved  using  a  finite  difference  discretization  of  the  time  interval 
[to,tf]  with  a  uniform  mesh  having  stepsize  At  and  points  fo,ti,  •  •  •  An  —  tf-  The  approximate  values  of  2  at 
these  times  are  denoted  by  Zo,  ■  ■  ■  ,Zn  .  A  central  difference  approximation  of  the  temporal  derivative  yields  the 
system 

^  izj+ 1  ~  zj ]  =  2  Zj)  +  F{tj. (_i,  Zj+i)\ 

E0z0  =  [2/o,  0]T  (40) 


EfZN  =  [0,  0]T 

for  j  =  0,  •  •  •  ,  N  —  1.  The  matrices  E0  and  Ef  represent  the  boundary  conditions  at  the  initial  and  final  times 
which  are  given  by  the  relations 
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Here  /  denotes  a  2 Nw  x  2NW  identity  matrix  where  Nw  denotes  the  number  of  basis  functions  employed 
in  the  spatial  approximation  of  the  state  variables  given  in  Section  3.2.  The  boundary  conditions  are  written 
in  this  form  to  ensure  the  initial  plate  displacement  and  velocity  are  fixed  according  to  the  prescribed  initial 
condition  in  (39),  while  the  solution  to  A (tf)  is  related  to  the  unknown  plate  displacement  and  velocity  at  tf 
which  must  be  determined  according  to  the  optimality  system. 

The  determination  of  a  solution  vector  Zh  =  [zo,  ■  ■  ■  ,  zjv]  hr  (40)  can  be  expressed  as  the  problem  of  finding 
Zh  which  solves 


H*h)  =  0 

where  the  system  of  equations  given  in  (42)  is  represented  by  T(zh)  €  R4(JV+1)JV“v 


(42) 


Hzh)  = 


F0 

Fi 


Fj 


Fn- i 
b(z0,zN) 


Fj  —  ^  [zj+ i  Zj\  ^  [F{tj,  Zj)  +  F(tj+i,  Zj+i)] 


b{z o,  zN)  =  E0z0  +  EfZN  - 


Vo 

0 


A  quasi-Newton  iteration  of  the  form  z^+1  =  zfc  +  where  solves 


(43) 


F\zkh)eh  =  ~F{zkh), 


(44) 


is  used  to  determine  the  solution  to  the  nonlinear  system  (42).  The  Jacobian  E'(z^)  has  the  structure 


F'(zh)  = 


So  R0 
Si  R[ 


Eo 
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Sn-i  Rn-i 
Ef 


where 


The  matrix  A(ti)  is  the  linearization 


which  yields  the  representation 


Si  ~  At1  2A^ 
Rl=  At1  ~  2A{fi+l) 


dF 

A(U )  =  —(ti,Zi) 


-  1 

I 

0 

1 

A 

£b[u*} 

Si  =  — — 

At 

0 

I 

_  2 

-Q 

-AT 

for  Si.  The  representation  for  Ri  is  similar. 


Direct  solution  of  (44)  is  infeasible  due  to  the  large  number  of  basis  functions  and  time  increments  required 
to  resolve  the  solution  over  a  reasonable  time  interval.  The  structure  of  the  Jacobian  can  be  employed  to  reduce 
both  memory  and  computational  requirements  to  the  level  of  solving  ANW  x  4 Nw  systems.  Following  [39],  this 
is  accomplished  by  expressing  the  Jacobian  in  terms  of  an  analytic  LU  decomposition 


T\zl)  =  LU 


(45) 


where 


So 


Si 


L  = 


E0  -E0(Sq1Ro) 


Sn-i  0 

N—2  N—l 

E0  (-1  )i(Sr1JRi)  Ef  +  E0  J]  (-lJ^Sr1^) 


i=0 


i= 0 
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The  direct  solution  of  the  lower  triangular  system  is  first  obtained  followed  by  direct  solution 

of  the  upper  triangular  system  U =  (£.  This  leads  to  the  solution  of  the  system  (44). 

Remark  1:  The  nonlinear  optimal  control  problem  includes  certain  approximations  to  improve  numerical 
efficiency.  For  example,  the  term  J^B[u*]  in  5)  and  Ri  is  not  trivial  and  is  assumed  to  be  linear  in  the  present 
model.  This  leads  to  a  suboptimal  nonlinear  controller  where  Si  =  S  and  J?j  =  R  for  all  time  steps.  Furthermore, 
the  condition  number  on  the  matrices  S  and  R  is  partially  dependent  on  the  time  step  used  in  the  simulations. 
In  the  simulations  presented  here,  the  time  step  At  =  0.01  s  was  employed.  Simulations  using  smaller  time 
steps  resulted  in  practically  the  same  dynamic  response  until  S  and  R  become  ill-conditioned  which  introduced 
a  numerical  instability.  This  occurred  when  At  <  0.005  s. 

The  term  B^(u*)  in  (26)  contains  nonlinearites  and  hysteresis  domain  switching  which  occurs  at  moderate 
to  high  field  levels.  This  value  can  be  computed  by  numerically  by  differientiating  the  constitutive  response 
during  the  control  simulation,  although  numerical  errors  may  introduce  instabilities.  Because  the  control  input 
is  considered  for  biased  minor  hysteresis  loops,  potential  numerical  instabilities  are  mitigated  by  assuming 
the  tangential  magnetic  susceptibility  is  equal  to  the  linear  coefficient  Xm-  Unmodeled  dynamics  resulting 
from  both  the  approximation  in  the  Jacobian  and  the  tangential  magnetic  susceptibility  are  addressed  through 
perturbation  feedback  control  presented  in  Section  4.3. 

In  Section  3,  a  reduction  in  the  number  of  control  inputs  was  implemented  in  the  control  design  by  assuming 
diametrically  opposed  magnetic  fields  on  each  top  and  bottom  actuator  pair  resulted  in  diametrically  opposed 
magnetization.  Although  this  simplifies  the  control  design,  nonlinear  H—M  behavior  may  give  rise  to  differences 
in  magnetization  under  moderate  to  high  magnetic  field  inputs.  In  the  present  model,  the  unmodeled  dynamics 
are  addressed  through  perturbation  feedback,  although  the  control  design  could  also  be  modified  to  include  the 
additional  degree  of  freedom  on  each  actuator  couple. 

Remark  2:  The  analytic  LU  decomposition  previously  developed  by  Smith  [39]  allows  good  spatial  and 
time  resolution  of  the  plate  dynamics.  When  considering  the  time  step  At  =  0.01  s  over  the  time  interval 
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[0.45,2.50]  s  and  spatial  approximation  with  4 Nw  =  100,  approximately  30,000  unknown  coefficients  must  be 
determined.  The  LU  decomposition  method  reduces  the  matrix  system  to  solving  100  unknowns  for  each  time 
step.  Numerical  stabilities  were  also  checked  for  larger  system  sizes.  The  time  step  was  decreased  to  0.005  s  and 
the  number  of  basis  functions  increased  to  8  (47VU,  =  324).  This  resulted  in  determining  over  130,000  coefficients 
which  was  possible  by  solving  the  324  coefficients  for  each  time  step  separately.  A  stable  solution  was  obtained 
that  was  practically  identical  to  the  one  where  41V,„  =  100  and  At  =  0.01  s  over  the  same  time  interval. 

4.3.1  Numerical  Example  -  Truncated  External  Force 

Enhancement  of  active  vibration  damping  is  significant  when  the  nonlinear  control  method  is  employed.  In 
Figure  7,  the  same  disturbance  load  (37)  used  in  the  linear  model  is  applied  to  the  Terfenol  plate  structure. 
When  to  =  0.45  s,  the  disturbance  load  is  again  set  to  zero  and  control  is  instantaneously  applied.  By  increasing 
the  penalty  on  the  states  through  Q  and  relaxing  the  penalty  on  u  through  R  as  given  in  Table  2  ,  the  control 
input  reaches  moderate  to  high  drive  levels.  Although  this  creates  nonlinearity  and  hysteresis  in  the  Terfenol 
transducers  similar  to  that  shown  in  Figure  5(b),  the  nonlinear  control  method  effectively  compensates  for  this 
behavior.  The  plate  vibration  at  the  point  [tx.  ty\  is  eliminated  in  less  than  half  the  time  relative  to  the  linear 
simulation  in  Figure  4. 


E 


E 


(a) 


(b) 


Figure  7:  Nonlinear  feedback  control  law  for  the  open  loop  control,  (a)  Plate  displacement  for  the  uncontrolled  ( - ) 

and  controlled  response  (^— ).  (b)  Relationship  between  magnetization  and  magnetic  field. 
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4.3.2  Numerical  Example  -  Periodic  External  Disturbances  and  Operating  Uncertainties 

The  effect  of  disturbances  on  the  plate  is  investigated  using  the  nonlinear  control  method  by  applying  a 
nonzero  G{t )  in  the  optimality  system  given  by  (39).  A  disturbance  load  of  25  N  is  again  applied  uniformly  to 
the  plate  at  the  previously  noted  frequencies.  Figure  8  illustrates  significantly  improved  vibration  attenuation 
when  the  nonlinear  method  is  employed.  Further  details  regarding  the  level  of  vibration  attenuation  attainable 
is  discussed  in  Remark  3  of  Section  4.4.1. 

The  previous  simulations  have  focused  on  nonlinear  open  loop  control  design  which  raises  the  question 
regarding  robustness  to  operating  uncertainty.  To  investigate  this  issue,  the  open  loop  control  trajectory  is 
applied  0.05  s  late.  In  addition,  the  multiple  frequencies  of  excitation  computed  in  the  open  loop  case  in 
(37)  are  perturbed  in  the  frequency  domain  by  +4  Hz.  In  the  simulation,  the  plate  is  initially  excited  at  the 
unperturbed  periodic  force  and  then  the  frequency  perturbation  is  applied.  Continuity  in  the  perturbed  and 
unperturbed  disturbance  forces  is  ensured  at  the  onset  of  turning  on  the  control  signal  through  the  disturbance 
force 


.  ,  [  25Ei=i  [sin(27r/f*)]  ,  t  <  .45 

g(t,x,y)=<  (46) 

[  25  Eli  [sin(27r/ft  -  1.8tt)]  ,  .45  <  t  <  2.5 

where  ff  are  the  expected  frequencies  given  in  Section  4.1.1  and  /“  are  the  ‘actual’  frequencies  shifted  by  +4  Hz. 
The  factor  1.87T  ensures  continuity  of  g. 

Figure  9  illustrates  complete  degradation  in  control  authority  when  operating  uncertainties  are  present. 
In  Figure  9(a)  the  effect  of  turning  on  the  control  signal  0.05  s  is  simulated.  It  is  illustrated  that  complete 
loss  of  control  authority  occurs  when  open  loop  control  is  applied.  In  Figure  9(b),  the  additional  effect  of 
perturbations  in  the  assumed  force  excitation  is  simulated.  In  this  case,  the  control  is  again  turned  on  0.05  s 
late  and  disturbance  loads  from  (46)  are  applied.  Complete  loss  of  control  authority  again  occurs  as  illustrated 
in  Figure  9(b). 
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Figure  8:  Performance  of  the  nonlinear  control  law  when  disturbance  loads  are  present,  (a)  Plate  displacement  for  the 
uncontrolled  ( - )  and  controlled  response  (b)  Relationship  between  magnetization  and  magnetic  field. 


Figure  9:  (a)  Performance  of  the  nonlinear  control  law  when  the  control  is  turned  on  0.05  sec  late,  (b)  The  additional 
effect  of  perturbed  periodic  disturbance  loads.  The  operating  uncertainty  causes  a  complete  loss  of  control  authority. 
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4.4  Perturbation  Control 


Open  loop  control  suffers  from  a  lack  of  robustness  to  operating  uncertainties  and  unmodeled  dynamics.  It 
has  been  shown  that  robustness  to  various  types  of  uncertainties  can  be  significantly  improved  by  considering 
perturbation  control  techniques  [5,  27].  The  feedback  control  input  is  determined  from  a  linearization  of  the 
system  about  the  optimal  control  pair  (u*(t),y*(t)).  Feedback  control  Su*(t)  is  designed  to  attenuate  pertur¬ 
bations  in  the  system  that  may  originate  from  external  forces  or  initial  conditions  described  in  the  previous 
section.  Since  the  theory  is  based  on  a  linearized  system  with  quadratic  constraints,  LQR  theory  can  be  utilized. 
This  provides  an  efficient  algorithm  for  implementing  the  control  design  in  real-time  by  first  computing  the  open 
loop  control  offline  and  using  online  feedback  of  Su*(t). 

The  modified  cost  functional  given  by  (22)  and  constraint  given  by  (20)  are  used  to  determine  the  pertur¬ 
bation  feedback  control.  Since  the  optimal  control  pair  (it*  (f) ,  y*  (t))  minimizes  the  first  order  variation  in  J, 
the  second-order  terms  are  expanded  to  determine  ( Su*(t),5y*(t )).  The  second-order  variation  in  J  and  the 
first-order  variation  in  the  constraint  are  given  by, 


and 
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(47) 


5y(t)  =  A6y(t)  +  BSu(t)  +  SG(t) 


(48) 


<ty(0)  =  Do 

where  Su  and  Sy  are  first-order  variations  about  u*  and  y* .  The  optimal  perturbation  control  Su*  is  that  which 
minimizes  (47)  subject  to  (48). 

For  the  Hamiltonian  given  by  (23),  the  second  variation  S2  J  is, 


82  J  =  ^  f  {8yTQ8y  +  5uTRSu}  dt  (49) 

2  ~'t0 

so  that  the  LQR  theory  outlined  in  Section  4.2  can  be  directly  employed  to  obtain  Su*(t)  and  Sy*(t).  The 
overall  control  for  the  system  is  then  taken  to  be  u*  ( t )  +  Su*  (t)  with  the  optimal  state  given  by  y*  (t)  +  Sy*  ( t ) . 


4.4.1  Numerical  Example 

The  robustness  of  the  perturbation  control  method  is  compared  to  the  previous  example  where  a  perturbed 
initial  value  and  disturbance  load  are  introduced.  As  illustrated  in  Figure  9,  perturbation  in  initial  conditions 
and  external  disturbances  is  sufficient  to  destroy  the  authority  of  the  open  loop  nonlinear  control.  The  two 
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previous  cases  are  again  considered.  First  perturbations  in  only  the  initial  value  are  considered  and  then  the 
additional  effects  of  perturbations  in  the  expected  periodic  disturbance  load  are  evaluated. 

When  perturbations  in  the  disturbance  load  are  absent,  the  control  law  is  given  by, 

Su*(t)  =  -R-'BFmyit)  (50) 

where  II  is  the  solution  of  the  algebraic  Ricatti  equation  for  the  perturbed  system. 

The  perturbation  in  disturbance  loads  is  represented  by  the  relation 

5 

Sg{t,  x,  y)  =  25  ^  [sin(27r/“<  -  1.8tt)  -  sin(27r/ft)]  (51) 

»= l 

.45  <  t  <  2.5  where  the  ‘actual’  frequencies  /“  have  again  been  shifted  by  +4  Hz  relative  to  the  expected 
frequencies  ff. 

The  control  law 


Su*(t)  =  -R~1Bt  [H8y(t)  -  Sr(t )] 


(52) 


compensates  for  the  disturbance  loads  where  Sr(t)  is  governed  by  the  following  auxiliary  equation, 


6r(t)  =  -[A-  BR~1BtU]T  r(t)  +  USG(t) 
5r(0)  =  8r(r) . 


(53) 


The  resulting  vibration  control  is  illustrated  in  Figure  10.  In  comparison  to  Figure  9(a),  the  effect  of 
turning  on  the  control  input  0.05  s  late  is  effectively  compensated  using  perturbation  feedback  as  illustrated  in 
Figure  10(a).  In  Figure  10(b),  the  additional  effect  of  perturbations  in  the  disturbance  load  is  simulated.  The 
perturbation  control  input  computed  using  (52)  is  shown  to  effectively  compensate  for  turning  on  the  control 
late  as  well  as  perturbation  in  the  disturbance  loads.  This  demonstrates  significant  enhancement  in  control 
authority  in  comparison  to  Figure  9(b).  It  is  noted  that  through  the  use  of  the  perturbation  feedback  control, 
vibration  attenuation  comparable  to  that  in  the  unperturbed  system  in  Figures  7  and  8  is  obtained. 

Remark  3:  The  open  loop  nonlinear  control  and  perturbation  feedback  simulations  given  in  Figures  8 
and  10(b)  do  not  completely  eliminate  plate  vibration.  This  is  not  a  deficiency  in  the  control  design,  but  a 
physical  limitation  of  the  Terfenol  actuators.  The  current  location  and  number  of  actuators  have  not  been 
optimized.  Improvements  can  be  made  by  determining  the  appropriate  location  and  number  of  actuators  for  a 
given  application. 
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Figure  10:  Performance  of  the  perturbation  feedback  control  law  when  operating  uncertainties  are  present, 
open  loop  control  is  turned  on  0.05  sec  late  and  zero  disturbance  loads  are  present, 
disturbance  loads  given  by  (51)  is  included. 


5.  Concluding  Remarks 

A  model-based  nonlinear  optimal  control  method  was  developed  to  actively  attenuate  plate  vibration  using 
magnetostrictive  materials.  While  linear  control  methods  provide  marginal  vibration  attenuation  when  the  input 
field  is  limited  to  the  linear  regime,  performance  improvements  can  be  obtained  by  driving  the  transducers  at 
higher  field  levels.  While  moderate  to  high  field  inputs  provide  the  potential  for  improved  vibration  control, 
the  resulting  nonlinear  and  hysteretic  constitutive  behavior  cannot  be  effectively  compensated  using  linear 
control  methods.  This  constitutive  behavior  was  directly  incorporated  into  the  optimal  control  design  which 
significantly  improved  vibration  control.  Moreover,  robustness  to  operating  uncertainties  was  obtained  by 
feedback  of  perturbations  around  the  optimal  trajectory,  thus  providing  a  hybrid  optimal  control  strategy  for 
real-time  applications. 

Whereas  the  nonlinear  control  design  for  vibration  attenuation  is  developed  in  the  context  of  magnetic 
transducers,  the  unified  nature  of  the  framework  permits  direct  extension  of  the  models  and  model-based 
control  methods  to  ferroelectric  and  ferroelastic  compounds  [40,  43,  44].  The  extension  of  the  control  theory 
to  specific  PZT  and  SMA  applications  and  experimental  implementation  of  the  control  designs  for  a  range  of 
ferrioc  transducers  operating  in  hysteretic  and  nonlinear  regimes  is  under  present  investigation. 

A  complementary  control  problem  is  that  of  tracking  a  reference  trajectory.  In  the  context  of  smart  material 
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actuators,  such  tracking  problems  include  high  accuracy  and  high  frequency  machining  of  out-of-round  ingots 
using  Terfenol-D  transducers  as  well  as  nanoscale  positioning  of  PZT  stages  in  an  atomic  force  microscope 
(AFM).  Whereas  the  form  of  complementary  ‘tracking  equations’  is  analogous  to  the  disturbance  differential 
equations  employed  here,  the  implementation  of  the  tracking  problems  differs  fundamentally  in  the  perturbation 
control  formulation.  The  development  of  a  complementary  nonlinear  optimal  control  design  for  high  accuracy, 
high  frequency  tracking  using  a  Terfenol-D  transducer  is  addressed  in  [30] . 
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Appendix 

We  summarize  here  the  moment-displacement  relations  used  in  formulating  the  thin  plate  structure;  details 
can  be  found  in  [40].  The  moment  components  are  given  in  terms  of  the  stiffness  and  damping  relations, 

int  Eph3  ( d 2w  d2w\  cph3  (  d2w  d3w  \ 

x  12(1  —  r-2)  p  dy2  )  12(1  —  is2)  \dx2dt  ^  pdy2dt) 

int  _  Eph3  ( d2w  d 2w\  Cph3  (  d3w  d3w  \ 
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Cph 3  /  d3w  \ 

24(1  +  Up)  \dxdydt J 


These  relations  are  substituted  into  the  weak  form  of  the  model  given  by  (15)  to  obtain  the  governing 
equations  in  terms  of  the  plate  displacement  w  and  the  test  functions  4>.  The  approximate  solution  to  this 
equation  is  obtained  by  substituting  the  plate  displacements  given  in  (17)  into  (15).  This  gives  rise  to  the  finite 
dimensional  formulation  in  (18).  The  resulting  mass,  damping,  and  stiffness  matrices  in  (18)  are  given  by 


Pp^j^kdu 


H-fc  =  \Ki]jk  +  [K2]jfc  +  [K3]jk  +  [K4]jfe  +  [K5]jk 


where 
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ijk 
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jk 


ijk 
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12(1  —  v2)/ 

/  dy2 
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VpEph 3 

\  d2$> 

d2<$>k 

12(1  —  v2)  / 

)  dy2 

dx2 

Eph 3  N 

\  d2*j 

d2$k 

dto 


dto 


dxdy  dxdy 


dto 


The  damping  matrix  C  is  constructed  in  a  manner  similar  to  the  stiffness  matrix  K.  The  disturbance  load 
is  defined  by 


[g ]j  =  J  O^jdoJ 

where  g  is  the  magnitude  of  the  disturbance  load  and  the  input  matrix  is 


£^cj)  •  ■ 

lhhi  =  /  I  +Xy(i)(x,y)-g~i- )  du. 


(54) 


in  \  nx-  dy2 

where  the  characteristic  functions  are  defined  in  (10)  and  (11).  This  ensures  consistency  with  the  applied 
moment  components  in  the  x  and  y  directions. 
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